Automatic estimation of point-spread-function 
for deconvoluting out-of-focus optical coherence 
tomographic images using information entropy- 
based approach 

Guozhong Liu, 1 ' 2 Siavash Yousefi, 1 Zhongwei Zhi 1 and Ruikang K. Wang 1 '* 

'Department of Bioengineering, University of Washington, 3720 15th Avenue Northeast, Seattle, Washington 98195, 

USA 

2 Department of Photo-electronic Information and Communication Engineering, Beijing Information Science and 
Technology University, Beijing 100192, China 
"wangrk@u. Washington, edu 

Abstract: This paper proposes an automatic point spread function (PSF) 
estimation method to de-blur out-of-focus optical coherence tomography 
(OCT) images. The method utilizes Richardson-Lucy deconvolution 
algorithm to deconvolve noisy defocused images with a family of Gaussian 
PSFs with different beam spot sizes. Then, the best beam spot size is 
automatically estimated based on the discontinuity of information entropy 
of recovered images. Therefore, it is not required a prior knowledge of the 
parameters or PSF of OCT system for de-convoluting image. The model 
does not account for the diffraction and the coherent scattering of light by 
the sample. A series of experiments are performed on digital phantoms, a 
custom-built phantom doped with microspheres, fresh onion as well as the 
human fingertip in vivo to show the performance of the proposed method. 
The method may also be useful in combining with other deconvolution 
algorithms for PSF estimation and image recovery. 
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1. Introduction 

Optical coherence tomography (OCT) is a non-invasive, micrometer resolution, three- 
dimensional optical imaging technique capable of capturing backscattering photons from 
within optical scattering media such as biological tissues [1, 2]. One of the advantages of 
OCT is that the axial and the lateral resolutions are decoupled from each other. The axial 
resolution of the OCT system is determined by the spectral bandwidth of light source and is 
defined as the half of the coherence length of the light source which can be expressed by 

R = 1-0.44^- (i) 

2 AA 

where l c is the coherence length, X 0 is the center wavelength and AX is the spectral full- 
width at half maximum (FWHM) of the light source. 

The lateral resolution of the OCT system is determined by the objective lens that is used to 
deliver/focus the light onto the sample and is expressed by 

=1.22 A ° (2) 

where NA obj is the numerical aperture (NA) of the objective lens which is inversely 
proportional to the lateral resolution. On the other hand, the depth-of-focus (DOF) is given by 

Z =2 ^- (3) 

DOF NA 2 

where n is the average refractive index of the sample. In general, the lateral resolution can be 
enhanced by the increase of numerical aperture. Unfortunately, the increase of the numerical 
aperture leads to the reduction of DOF. The relationship between the lateral resolution and 
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DOF is schematically illustrated in Fig. 1 for the cases of low (left figure) and high NAs (right 
figure). From the illustrations, we can appreciate that only the portion of the 3D image within 
DOF exhibits the desired lateral resolution, while that in the de-focal plane is blurred. There 
exist several dynamic focusing methods in the literature to mitigate the compromise between 
lateral resolution and DOF. In general, these methods can be divided into hardware-based and 
software-based (digital) approaches. 



Low Numerical Aperture High Numerical Aperture 

x x 




(a) (b) 

Fig. 1. Schematic illustration on the effect of numerical aperture on the desired lateral 
resolution of OCT images, in the case of (a) low NA, and (b) high NA. 



1.1 Hardware-based approaches 

Lexer et al. [3] investigated an optical setup to shift the focus of the beam illuminating the 
object through the object depth without changing the path length in the corresponding 
interferometer arm and achieved lOrim lateral resolution over an object depth of 430 urn in 
human cornea. Qi et al. [4] designed a high-speed dynamic focus control system based on a 
micro-electro-mechanical mirror, in which the silicon nitride deformable mirror shifted the 
focus of the sample beam of the OCT interferometer in synchrony with the coherence-gate 
scan in the reference arm. As a result, the coherence gate remained at the beam focus during 
the whole imaging process. Divetia et al. [5] used a 1 mm liquid-filled polymer lens for 
endoscopic OCT applications that dynamically provided scanning depth of focus by changing 
the hydraulic pressure within the lens which enabled dynamic control of the focal depth 
without the need for articulated parts. This configuration was shown to have resolving power 
of 5 urn over an object depth of 2.5 mm. Xie et al. [6] developed a gradient index (GRIN) lens 
rod based probe for endoscopic spectral domain OCT with dynamic focus tracking. The focus 
of their endoscopic OCT probe could be dynamically adjusted at 500 mm/s without changing 
reference arm length to obtain high-quality OCT images for contact or non-contact tissue 
applications. The dynamic focusing range of the probe was achieved between 0 to 7.5 mm 
without moving the probe beam itself, while the imaging depth was 2.8 mm. To overcome the 
trade-off between lateral resolution and focusing depth, Ding et al. [7] incorporated an axicon 
lens with a top angle of 160° into the sample arm of an interferometer and maintained lOum 
lateral resolution over a focusing depth of 6 mm. Holmes et al. [8] presented a multi-beam 
OCT system that overcomes the problem of limited lateral resolution inherent in single-beam 
Fourier domain OCT and reduces speckle noise. However, the hardware-based methods 
require special hardware set-up configuration in the system design which can limit the 
scanning speed and raise the cost of the OCT system. 

1.2 Digital approaches 

There also exist several digital methods in the literature to compensate image degradation out 
of the focal depth zone by solving inverse scattering and deconvolution problems. Ralston et 
al. [9] proposed a novel computational image-formation technique called interferometric 
synthetic aperture microscopy (ISAM) by solving the inverse scattering problem for 
interference microscopy. ISAM could achieve reconstructed volumes with a resolution in all 
planes that was equivalent to the resolution achieved only at the focal plane for the 
conventional high-resolution microscopy. Yu et al. [10] developed a two-dimensional scalar 
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diffraction model to simulate the wave propagation process from out-of-focus scatters within 
the short coherence gate of the OCT system and high-resolution details could be recovered 
from outside the depth-of-field region with minimum loss of lateral resolution. These two 
methods account for the diffraction and the coherent scattering of light by the sample, thus are 
the most accurate representation of the sample's microstructure outside the focus region, with 
ISAM being the most elegant one. However, they are extremely sensitive to the phase 
stability of the measurements, particularly when the in vivo imaging is involved where the 
sample movement is inevitable. Therefore, there is still a need to seek for alternatives to 
mitigate this issue. In this pursuit, the best solution is to use the amplitude information of the 
OCT signals because it is insensitive to the optical phases of the system. Over the last decade, 
efforts have been paid that explore the amplitude information of the OCT signals for the 
purpose of de-blurring the OCT images. 

Yasuno et al. [11] applied a numerical decon volution method to cancel lateral defocus in 
Fourier-domain OCT using a depth-dependent lateral PSF and applied the method to OCT 
images of in vivo human anterior eye segment. Kulkarni et al. [12] proposed a linear shift- 
invariant system model to describe coherent light specimen interactions in OCT images and 
applied an iterative deconvolution algorithm to enhance the sharpness of the images of 
biological structures. A constrained iterative restoration algorithm was utilized to estimate the 
impulse response of the system model using a prior knowledge of the properties of that 
response [13]. The algorithm was applied to a defocused OCT image of fresh onion sample 
and the deconvolved OCT image revealed polygonal cellular structure in the specimen. 
Ralston et al. [14] proposed a method for reducing transverse blurring and improving the 
transverse resolution in OCT using Gaussian beam deconvolution. First, the direct inverse 
problem was investigated using two types of regularization, truncated singular value 
decomposition and Tikhonov. Then, an iterative expectation-maximization algorithm, the 
Richardson-Lucy algorithm, with a beam-width-dependent iteration scheme was developed. 
Using a dynamically iterative Richardson-Lucy algorithm reduced the transverse blurring by 
providing an improvement in the transverse PSF for sparse scattering samples in regions up to 
two times larger than the confocal region of the lens. Also, Wiener and Richardson-Lucy 
OCT image restoration algorithms were implemented and compared with each other where 
the images deconvolved using Richardson-Lucy had a better contrast and quality [15]. 
Besides, JM Schmitt [16] derived CLEAN deconvolution kernel to provide an effective way 
of improving resolution and contrast in optical coherence tomography. Rolland et al. [17] 
developed a Gabor-based fusion technique based on the concept of the inverse local Fourier 
transform and the Gabor's signal expansion to produce a high lateral resolution image 
throughout the depth of imaging. 

Since these previous digital methods require a prior knowledge of the optical parameters 
or PSF of the optical system, an initial calibration experiment is required. Tomlins et al. [18] 
and [19] made PSF phantoms doped with a small quantity of calibrated spherical micro- 
particles to estimate the Gaussian PSF parameters and researched the potential PSF degrading 
effects of optical dispersion and spherical aberration. Tomlins et al. [20] also adopted 
femtosecond laser subsurface micro-inscription techniques to fabricate an OCT test artifact 
for validating the resolution performance of a commercial OCT system. The method [18] and 
[19] were generalized by Woolliams et al. [21] to estimate 2D PSF at multiple locations 
across the B-scans to account for the variation in PSF as a function of both depth and lateral 
position. Agrawal et al. [22] present a PSF measurement approach using a nanoparticle- 
embedded phantom, towards the development of standardized test methods. However, the 
optical parameters or the PSF of OCT system needed to be calculated by performing a 
separate experiment on a specialized phantom to calibrate the parameters and the parameters 
may not be valid for complex heterogeneous structures. 

In this paper, by using the forward OCT model previously developed by Ralston et al [14], 
we propose an automatic PSF estimation method based on discontinuity of information 
entropy and Richardson-Lucy deconvolution algorithm to improve the image contrast and 
quality. The advantage of this method is that there is no need to implement a new optical 
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hardware or to measure the optical parameters of the OCT system, and importantly it is 
relatively in-sensitive to the phase-instability of the measurements. Since the Richardson- 
Lucy algorithm tends to concentrate energies near boundaries, it provides a good 
approximation to cellular boundaries and subcellular features and is more robust to noise. 
Therefore, combining PSF estimation method and Richardson-Lucy algorithm can estimate 
PSF of OCT at different defocal planes automatically and objectively while decreasing the 
transverse blurring. 

2. Image recovering principle and algorithms 

In Gaussian optics model, the lateral PSF of OCT is approximated by [14]: 

2 , 2 

h(x,y) = exp(-2—^f-) 
where W z is the depth-dependent beam spot size given by 

W z =W 0 [l + (— ff 

Zr ^ 

An 

where W 0 is the minimum value of W z which occurs at z = 0, z is the axial coordinate and 

z = z R at the boundary of the confocal region (Rayleigh range) . The assumption in (5) is that 

Gaussian beam propagates in a direction normal to the lens and the refractive index n is 
constant along the beam. However, the beam profile deviates from the ideal model and the 
transverse resolution is not constant because index of refraction is not constant in the sample. 
Also, the beam profile depends on the depth and the focal region of the lens. 

In OCT system, the final image at each depth location I(x, y) is the convolution of the 

lateral PSF h(x, y) and the sample profile 0(x, y) which can be formulated by 

l(x,y) = 0(x,y)* h(x, y) + n = jJ^O(x-x 0 ,y- y 0 )h(x 0 ,y 0 )dx 0 dy 0 +n (7) 

where n is additive noise and tissue heterogeneous noise [23]. By estimating the best beam 
spot size W z at each depth location, the lateral resolution at out-of-focus depth locations can 
be improved by solving an inverse problem. 

2.1 Image recovering flow diagram 

The flow diagram of the proposed method for recovering defocused FD-OCT images is given 
in Fig. 2. 

Stepl: The non-linear spectral interferogram captured in an A-scan is rescaled into linear 
k-space followed by Fast Fourier Transform (FFT) to become A-line complex information 
(phase and amplitude) along the z axis, whose amplitude represents the structural information. 
The 3D structure of a sample is acquired by two-dimensional scanning of the probe beam by a 
paired galvanometer mirror. 

Step2 and Step3: The 3D structure image is resampled along the depth (in the z direction) 
to obtain a sequence of en-face images I(x,y) at different depth locations followed by a 
Gaussian low-pass filter. Here, the Gaussian LPF acts to reduce the noise levels in the image, 
including the system noise and the spurious effects from speckle, so that the recovered image 
becomes smooth. It should be noted that applying a Gaussian LPF will reduce the lateral PSF 
of OCT system as follows: 



(4) 

(5) 
(6) 



#150751 - $15.00 USD Received 8 Jul 201 1; revised 14 Aug 201 1; accepted 17 Aug 201 1; published 31 Aug 201 1 
(C) 201 1 OSA 12 September 201 1 / Vol. 19, No. 19 / OPTICS EXPRESS 18139 



I' = f*I = f*(h*0 + n) = f*h*0 + f*n«(f*h)*0 = h'*0 



(8) 



where / is the impulse response of Gaussian LPF, / is the raw image, h is the lateral PSF, O 

is the object being imaged (sample), n is additive noise, /' is the filtered image and h' is the 
lateral PSF after low-pass filtering. In the experiments, the low-pass filter was empirically 
selected such that its bandwidth is approximately 5 times broader than that of the system PSF. 



Step 1 



3D Spectral interferogram 
1 



Rescalina +FFT 



Recovering the original image 
using Richardson-Lucy algorithm 
and estimated PSF ( W z ) 



Step 2 Resampling in the z direction 



I 



Step 3 Gaussian low-pass filter 



PSF ( W.) estimation using 
Laplacian of evaluation function 

f • 



Richardson-Lucy deconvolution 
using a family of W z s 



Step 6 



Step 5 



Step 4 



Fig. 2. Flow diagram for recovering defocused OCT images. 

Step4 and Step5: Richardson-Lucy deconvolution method requires a prior knowledge 
about the lateral PSF, h(x, y) , which is a function of the actual beam spot size W . In order 
to estimate the best choice of W z at each depth location, the filtered en-face images I'(x,y) 
are deconvolved using Richardson-Lucy algorithm by a Gaussian family PSFs with stepped 
beam spot size ff, e[ff d ,ff l2 ] and an evaluation function defined in subsection 2.3 is 
analyzed. The evaluation functions tends to have a discontinuity around W - W m , so, 
PSF(The best choice of W ) is acquired by searching a peak of the Laplacian of evaluation 
function . 

Richardson-Lucy is an iterative method based on the maximum-likelihood estimation for 
Poisson statistical data, which is updated by 



O n+l {x,y) = O n (x,y)[h\-x-y)*- 



I'(x,y) 



h (x, y)*O n {x,y) 



(9) 



where I'(x, y) is the degraded image after Gaussian LPF, O n (x, y) is the estimate of the 
original image after n iterations, h'(x, y) is the modified lateral PSF and h'(x, y) * O n (x, y) is 
referred to as the reblurred image. The iterations will stop when the difference between O n+1 
and O n goes below a threshold. 

Step6: Each original en-face image can be recovered using Richardson-Lucy algorithm 
and estimated lateral PSF ( W z ). 

2.2 Discontinuity of recovered image around W z = W m 

The output of OCT imaging system is described mathematically by convolving its PSF with 
an input image. Let's assume that a sufficiently small object can be modeled as a Dirac delta 
function. Then, the output of the imaging system is its lateral PSF as follows: 

x 2 + 2 

I(x, y) = 0(x, y) * h(x, y) = 5(x, y) * h(x, y) = h(x, y) = exp(-2- — f-) (10) 



where / is the received image, O is the object profile, h is the lateral PSF and W za is the 
actual beam spot size. When received image / is deconvolved using different beam spot size 
W, , it can be described as follows: 
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2 2 

I(x, y) = 0'(x, y) * h'(x, y) = 0'(x, y) * exp(-2^-^) (11) 

W z 

where 0'(x, y) is the recovered image when beam spot size is W z . Equation (12) is obtained 
from Eq. (10) and Eq. (11) and 0'(x, y) can be solved by using 2-dimentional convolution 
theorem as Eqs. (12) and (13). 

2 2 2 2 

exp(-2^^) = C>'(x,y)*exp(-2 X + y ) (12) 
W W 



0'(x,y) = 2 J exp(-2 _ 2 * rl ) (13) 



W[ x 2 + y 2 

; ^eXp(-2 - 

W UW -W ) W -W 

zct V z za / z zci 

It can be observed that the recovered image 0'(x, y) will have a discontinuity at W z = W za . 
According to sampling property of delta function, any input image structure can be 
modeled as weighted superposition of delta functions. 



°(*. y) = SL Sli - n, y - m) (14) 



where A nm is the weight (intensity) of the image at each n and m location. Based on the 
definition of linearity, the results from a delta function can be generalized to this case. 
Therefore, the recovered image 0'(x, y) is also discontinuous at the point W z = W za - This is 

exactly the theoretical basis to estimate the actual W za . In order to find the discontinuity of the 
recovered image around W z = W za , an evaluation function based on information entropy is 
defined and the discontinuity of this evaluation function is used to estimate the actual W za . 

2.3 Evaluation function for PSF estimation 

Information entropy of the recovered image 0'(x, y) and the blurry out-of-focus raw image 
I(x, y) can be expressed by 

= -^ a P(0'(a)).log(p(0'(a))) (15) 
E W = -E* p(I(b)).\og(p(I(b))) (16) 

where p(0'(a)) is the marginal probability of the image O' , which is the probability that the 

intensity value of the image is equal to a. The joint entropy of O' and / images, which is an 
measure of uncertainty associated with them, can be expressed by 

E(0',I) joint =-Y Jab p{IO\a,b)).\og(p(IO\a,b))) (17) 

where p(IO'(a,b)) is the joint probability distribution of O' and / images, which is the 
probability that the intensity value of the image O' is equal to a and the intensity value of the 
image / is equal to b, respectively. Also, the mutual information entropy of the O' and / 
images can be expressed by 

E(0',I) mulual =E(0') + E(I)-E(0',I) Mnt (18) 
The evaluation function is defined by 

F(W 7 ) = (19) 

^ E{0\I) hml -E(0') 
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which is a function for W e [W.j,W, 2 ] . Because of the discontinuity of the recovered O' , the 

evaluation function tends to have discontinuity at W z = W za . The discontinuity of the 
evaluation function is equivalent to a peak in the Laplacian (second derivative) of that 
function which is used to locate the discontinuity of the evaluation function. Typical curves of 
evaluation function and Laplacian of evaluation function are displayed in Fig. 3. 




Fig. 3. The peak of the Laplacian of evaluation function (b) corresponding to the discontinuous 
point of evaluation function (a) is used to locate the actual beam spot size W . 



2.4 Image recovery 

If the object en-face image 0(x, y) is located in the focal plane (W z = W 0 ), the received image 
will be given by 



x 2 + 2 

I 0 (x, y) = 0(x, y)*h 0 (x, y) = 0(x, y) * exp(-2 * / ) 



(20) 



Assuming that the object en-face image 0(x,y) is located in the out-of-focus plane (W z 
W^, the received image will be given by 



x 2 + 2 

I(x, y) = 0(x, y) * h(x, y) = 0(x, y) * exp(-2 ) 



(21) 



The focused image Io(x,y) can be recovered from out-of-focus image l(x,y) by filtering the 
out-of-focus image with the compensation function g(x,y) as follows 

A) (*> y) = 8(x, y) * I(x, y) = g(x, y) * 0(x, y) * h(x, y) = 0(x, y) * h 0 (x, y) (22) 

which yields a compensation function for image recovery 



>(x, y) = 2 



W; 



w 0 2 (K 



■w 2 ) 



2 , 2 

exp(-2 — 5 — J —) 
W 2 -W 2 



(23) 



3. Simulation studies 



In order to show the performance of the proposed method, we performed a series of 
simulations on digital phantoms. First, a Dirac delta function digital phantom was generated 
and manually defocused at different W za . For each W za scenario, Richardson-Lucy algorithm 
was performed to recover the image on a family of W z values and the discontinuity (peak of 
the Laplacian) of the evaluation function was measured to estimate the beam spot size. Figure 
4(a) shows a 3D plot of the evaluation function for different W za s where the x axis is W z , the y 
axis is W za and z axis (intensity) is the evaluation function. It can be observed that the 
discontinuity of the evaluation function is at W z = W m . Also, the Laplacian of the evaluation 
function is shown in Fig. 4(b) where x axis is W z , y axis is W za and z axis is the Laplacian of 
the evaluation function. The peak values are observed around W z = W za , which was expected. 
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Fig. 4. When digital phantom was defocused at different W m , (a) the discontinuous point of 
evaluation function approximately equals W m for all W m s, and (b) the peak of Laplacian of the 
evaluation function approximately equals W m for all W m s. 




Fig. 5. The performance of the proposed method on a digital phantom, (a) Original image, (b) 
Out of focus image. Out of focus image was then degraded with (c) speckle noise, (d) Poisson 
noise and (e) both speckle and Poisson noise. Recovered images: (f) without the presence of 
noise. With the presence of different noise conditions, (g) speckle noise, (h) Poisson noise and 
(i) both speckle and Poisson noise. 

Figure 5(a) shows another digital phantom to simulate the cellular structures of onion. The 
image was manually defocused at W = 3(Fig. 5(b)) and the proposed method was performed 

to recover the image where the recovered image is shown in Fig. 5(f). In order to study the 
effectiveness of the proposed method in suppressing noise, the image was degraded by 
Poisson and speckle noise. Figure 5(c-e) show defocused and degraded images by speckle 
noise, Poisson noise, and speckle and Poisson noise, respectively. Figure 5(g-i) show the 
recovered images from different noise conditions. The multiplicative speckle noise was a 
zero-mean uniformly distributed random noise with variance 0.05. 

The experiment was repeated for different defocused depth locations ( W m ) and the 
difference between the actual and estimated a number of beam spot size (estimation error) 
was measured at different noise conditions. Figure 6 shows the error in estimating actual W 

under different noise conditions. It can be observed that the algorithm can effectively suppress 
the noise effect and the results are comparable with the experiment when no noise is present. 
The proposed method was specifically more effective when the image was slightly defocused 
(W g [4 9]). In general, the estimation error was acceptable and the proposed method could 
effectively estimate the beam spot size at the present of different noise conditions. 
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Fig. 6. Estimation errors when digital phantom was manually defocused at different W M for 
different noise conditions. The maximal absolute error is about 0.2 within the range 1 tol2 of 
W m . 

4. Experimental results 

The OCT system used to achieve 3D volume data is shown in Fig. 7, which is similar to one 
of previous studies [24]. Briefly, the system used a super luminescent diode light source, with 
a central wavelength of 1300 nm and a bandwidth of 80 nm that provided a 12 um axial 
resolution in air. The light is split into two paths in a fiber based Michelson interferometer. 
One beam is coupled onto a stationary reference mirror and the second beam is focused with a 
collimating lens and a focal lens (NA = 0.05) at the sample with a theoretical lateral 
resolution of -16 (im and a depth of focus of -350 urn. The focal spot on the sample was 
scanned using a pair of galvanometer mirror mounted in the sample arm. The spectral 
interferogram between the reference light and the light backscattered from the sample was 
sent to a homebuilt spectrometer via an optical circulator. The spectrometer consisted of a 
transmission grating (1200 lines/mm), a camera lens with a focal length of 100 mm, and 1024 
element line scan InGaAs detector (capable of a line rate of 47 kHz), with a wavelength 
coverage of 200 nm. A pilot laser is also coupled into the interferometer for beam guidance. 

In order to show the performance of the proposed method in real OCT images, three sets 
of experiments were performed on: 1) a phantom doped with microspheres that simulate the 
PSF of the system, 2) a fresh onion to show the efficiency of the proposed method on the still 
sample, and 3) in vivo tissue sample to demonstrate whether the proposed method has the 
potential for in vivo imaging applications. For the results presented below, the processing 
time of one frame image (512*512) was 1.23s for a HP dv6 laptop computer (processor: Intel 
i3 2.4GHz; RAM: 6G) when number of iterations was 10. Also note that during the image 
reconstruction process, there was no threshold applied to the images. Only when the OCT 
images were displayed, a fixed threshold (equal to the noise level in the raw OCT image) was 
applied to better show the final images. 
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Fig. 7. Schematic of the OCT system used in this study 

First, a phantom was made from gelatin, doped with microspheres, having a predominant 
diameter of 300 nm. Since the microspheres are small relative to the system resolutions, they 
can be approximated by delta functions and the phantom can be modeled as a superposition of 
delta functions. 3D OCT volume imaging was obtained from the phantom and the focus of the 
system was set to the center of the phantom. 

Figure 8(a) shows one B-scan of the phantom where the center of the image is in focus but 
the upper and lower portion of the image has a degraded lateral resolution and the 
microspheres look blurry. Then, the proposed automatic PSF estimation and image recovery 
(see Eq. (22) and (23)) method was applied on to the 3D volume data at different depth 
locations and the recovered image of the same B-scan is shown in Fig. 8(b). Note that the 
blurred image in the out-of-focus area is now focused after performing the proposed algorithm 
on the experimental data. In addition, we can also recover the raw B-scan image (Fig. 8(a)) to 

object profile (Fig. 8(c)), which are delta functions when microspheres are sufficiently small, 
by deconvolution (see Eq. (7)) using estimated beam spot size W m at different depth locations. 

Note that some connected microspheres in raw B-scan image are seperated and the image 
quality is further improved. 
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Fig. 8. (a) Original image where the focus is in the center and the upper and lower part of the 
image is out of focus, (b) De-blurred image, (c) Recovered image knowing that each point is a 
delta function. 



Then, an experiment was performed on a fresh onion and 3D OCT images with 256 A- 
scans and 256 B scans were acquired from the focus area and 1.5 mm away from the focal 
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point. Figure 9(a) shows one en face image of the onion in the focal depth and the proposed 
method is applied to improve the contrast and resolution of the image in the presence of noise. 
Figure 9(b) shows the recovered image from Fig. 9(a). Then, while the probe beam was kept 
unchanged, the sample was shifted downwards by 1.5 mm. Figure 9(c) is the resulted image 
that corresponded to the same depth location of sample shown in Fig. 9(a). The proposed 
method was applied to recover that image (Fig. 9(c)), and the result is shown in Fig. 9(d). 
Note the similarity between Fig. 9(d) and Fig. 9(b) shows onion structures at the same depth 
location. 




Fig. 9. (a) The en face image of the onion layer in the focal depth, (b) the defocused same en 
face image with (a) when the sample was shifted downwards by 1,5 mm, (b) and (d) are 
recovered noise-suppressed images from (a) and (c) respectively. The axial and lateral 
resolution is 12 pm and 16 um respectively. 



The recovered images of onion layers at different depths are shown in Fig. 10 where Figs. 
10(a-d) show the defocused onion layers at 1500 um, 1604 um, 1725 um and 1823 um, 
respectively. Figures 10(e-h) show the recovered images at their corresponding depth 
locations in the onion. 




Fig. 10. (a-d) Out-of-focus noising images at different depth locations: 1500 um, 1604 um, 
1725 um and 1823 um respectively, (e-h) Recovered images from the corresponding depth 
locations. The axial and lateral resolution is 12 um and 16 um respectively. 



To show the discontinuity during the recovering process, Fig. 1 1 gives the 2D plots of PSF 
estimation curve for recovering from Fig. 10(a) to Fig. 10(e), where the changes of evaluation 
function and Laplacian of the evaluation function with the beam spot size W z are displayed. It 
can be observed that the discontinuity of the evaluation function and peak value of the 
Laplacian of the evaluation function occur when W z = 4.0, meaning that the estimated beam 
size was 4.0 for this en face image. 

The theoretical beam spot size W, a is compared with the estimated beam size W, a in the 
onion at different depths and is shown in Fig. 12 where the x axis is the depth in um and y 



#150751 -$15.00 USD Received 8 Jul 201 1; revised 14 Aug 201 1; accepted 1 7 Aug 20 1 1 ; published 3 1 Aug 2011 
(C) 201 1 OSA 12 September 201 1 / Vol. 19, No. 19 / OPTICS EXPRESS 18146 



axis is the beam spot size in um. It can be observed that the estimated beam spot sizes agree 
with the expected theoretical values. The maximum estimation error is 4.2 um and the 
maximum relative error is 12.8% for a depth location between 1.5 mm to 2.3 mm. 



FCUfe) F '( Wz , 




Fig. 11. PSF estimation curve for recovering the image from Fig. 10(a) to Fig. 10(e): (a) the 
change of the evaluation function with the beam spot size W z and its discontinuity (pointed by 
arrow) and (b) the corresponding Laplacian of the evaluation function with the beam spot size 
W z and its peak value. 
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Fig. 12. The theoretical beam spot size at different depths (dashed line) in the fresh onion 
sample is compared with estimated beam spot size (solid line). The maximum estimation error 
is 4.2 pm and the maximum relative error is 12.8% within a depth from 1.5 mm to 2.3 mm. 

Lastly, the preliminary experiments were performed on the OCT images that were 
acquired from the human fingertip using the system described in Fig. 7 to demonstrate 
whether the proposed method is applicable in vivo. For this set of experiments, the system 
was running at 120 frames per second. 3D OCT images consisted of 256 A-scans and 200 B 
scans, covering an area of 1.5x1.2 mm 2 on the skin surface. The data was first obtained from 
the focus region. And then with the other system setup fixed, the finger was moved 500 um 
away from the focal point, and another 3D image was acquired. Because of the system spatial 
resolution was about 16 itm, the image quality enhancement for the regions within dermis was 
not easily appreciated. For this reason, we decided to show here the ability of proposed 
method to provide enhanced imaging quality for the de-focused sweat glands. Figure 13(a) 
shows one typical en-face image located at the focus region, which was dissected from the 3D 
volume data set. This en-face image represented a plane approximately cut through the 
upward oriented sweat glands with 90 degrees, in which the sweat glands (shown as the bright 
spots) and the skin ridges are clearly seen. Figure 13(b) shows one defocused en-face image at 
the approximately the same location as in Fig. 13(a). The proposed method was then used to 
de-blur the de-focused image. The result is given in Fig. 13(c), where it can be seen that the 
recovered image has qualitatively the same quality as the focused image. Quantitatively, the 
estimated entropies for the three images were 6.9699 (Fig. 13(a)), 7.0128 (Fig. 13(b)), and 
6.9356 (Fig. 13(c)), respectively. Thus, these preliminary results demonstrate that the 
proposed approach is efficient in de-blurring the in vivo OCT images. Further refinement of 
the current model for in vivo imaging applications are still in progress. 
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Fig. 13. The en-face images obtained from fingertip of a human volunteer, showing a plane cut 
through the sweat glands at -90°. (a) Image in focal area, (b) Defocused image after moving 
-0.5 mm away from the focal point, (c) Recovered image from (b). 

5. Conclusion 

We proposed an automatic PSF estimation method for recovering defocused OCT images 
based on the discontinuity of information entropy. A family of beam spot sizes were used to 
recover the original image using Richardson-Lucy algorithm and the maxima in the Laplacian 
of the introduced evaluation function was utilized to find the best PSF. The performance of 
the proposed method was first tested on digital phantoms and then on a custom-built phantom, 
onion cells as well as on human finger in vivo. The quantitative performance of the method 
was quantized by knowing the parameters of the system and comparing the estimated 
parameters with their real values. It was observed that the proposed method is effective to 
recover defocused images in the presence of noise and extend depth of imaging range to 
2.8mm. This method can be used in recovering defocused OCT images where the PSF is not 
known, especially when the index of refraction at different depths of the sample is unknown 
and not fixed. The proposed method used an approximated OCT forward model when NA of 
OCT system is relatively small. In future, the combination of the proposed PSF estimation 
method with accurate OCT forward model, such as that described in interferometric synthetic 
aperture microscopy [9], will be explored to account for the diffraction and the coherent 
scattering of light by the sample, so that a better representation of the sample's microstructure 
outside the focus region is presented. The proposed image recovery approach assumes that the 
PSF exhibits linear shift invariance throughout an en face plane, which was achieved by 
confining the beam scanning to the central region of the objective lens. 

The proposed PSF estimation method can be combined with other deconvolution 
algorithms such as wiener filters for PSF estimation and image recovery. Furthermore, the 
proposed PSF estimation and image deconvolution method is not limited to OCT and may be 
useful for image defocusing in other confocal optical system. 
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